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The equation of state of QCD matter for temperatures near and above the quark-hadron transition 
(~ 165 MeV) is inferred within a Bayesian framework through the comparison of data from the 
Relativistic Heavy Ion Collider and from the Large Hadron Collider to theoretical models. State-of- 
the-art statistical techniques are applied to simultaneously analyze multiple classes of observables 
while varying 14 independent model parameters. The resulting posterior distribution over possible 
equations of state is consistent with results from lattice gauge theory. 


I. INTRODUCTION 

Relativistic heavy ion collisions have been proposed as 
a means for investigating the equation of state of hot 
matter. For fixed target energies of E/A < 10 GeV, 
analyses of heavy ion collisions have significantly con¬ 
strained the compressibility of dense hadronic matter |T] 
for temperatures < 100 MeV. Higher energy collisions 
probe conditions near and above the transition temper¬ 
ature, where lattice calculations have shown that in a 
narrow temperature band, 150 < T < 200 MeV, the 
scalar quark condensate melts j2j, the degrees of freedom 
change [3], and the speed of sound has a minimum [4]. 
In fact, for some time the transition was postulated to 
contain a first-order phase transition accompanied by a 
sizable latent heat. 

In contrast to the progress of lattice calculations, ex¬ 
perimental determination of the equation of state at 
high temperature has remained semi-quantitative. The 
stunted progress has not been due to a shortage of ex¬ 
perimental observables that are known to be sensitive 
to the equation of state. Van Hove associated the de¬ 
pendence of the mean transverse momentum, (p t ), as a 
function of multiplicity as tool for determining the equa¬ 
tion of state [5p Two-particle femtoscopic correlations 
were proposed as a signal for a first-order phase transi¬ 
tion [B], Measurements of azimuthal elliptic flow, which 
are now mainly associated with determining the viscosity, 
were also shown to be sensitive to the equation of state 
Hi- Multiplicities, which are related to entropy, have 
also been used to constrain the equation of state [5j. Al¬ 
though femtoscopic analyses have shown that a first order 
equation of state with a large latent heat is highly un¬ 
likely (ID], and that an extremely stiff equation of state, 
such as that of a pion gas, is also inconsistent with data 
m, a more quantitative statement of how well the equa¬ 
tion of state is constrained has proven elusive. Even if 
analysis of experimental data cannot compete with lat¬ 
tice calculations in determining the equation of state for 
perfectly equilibrated matter, constraining the equation 
of state by experiment can help validate the statement 
that the matter created in heavy-ion collisions behaves 
like an equilibrated quark gluon plasma. 


The road block to turning these sensitivities into a 
more robust and rigorous determination of the equation 
of state has been the intertwined dependencies between 
the many unknown features and parameters of the model, 
and the numerous classes of measurement. Two develop¬ 
ments now make this next step possible. First, the mod¬ 
els used to describe the bulk behavior have converged to 
a standard framework based on relativistic viscous hy¬ 
drodynamics for the evolution of the high temperature 
region, > 165 MeV, m coupled to a microscopic simu¬ 
lation of the lower temperature hadronic stage based on 
binary collisions. The initial evolution, which feeds into 
the hydrodynamic description, remains rather undefined, 
but one can represent those uncertainties parametrically. 
The second development is in the statistical methodolo¬ 
gies and tools required to compare heterogenous data to 
models where a large number of parameters are required 
to encapsulate the many model uncertainties. Here we 
use the statistical tools described in |T3] to constrain 14 
parameters via a Markov-chain Monte Carlo. The statis¬ 
tical tools are based on a Gaussian-process model emula¬ 
tor, which allows one to estimate observables for a given 
point in parameter space by interpolating from a fixed 
number of full-model runs. 


II. METHODOLOGY 

Here we report on comparisons of model calculations to 
data from Au+Au collisions from the highest RHIC (Rel¬ 
ativistic Heavy Ion Collider) energy, 100A GeV + 100A 
GeV, and from Pb+Pb collisions at the Large Hadron 
Collider (LHC), 1.38A TeV + 1.38A TeV. The hydrody¬ 
namic and hadronic simulations were the same as those 
used in m to analyze RHIC data. The analysis in¬ 
volves 14 parameters, two of which vary the equation of 
state. The statistical method returns a sampling of the 
14-dimensional space that is weighted by the likelihood, 

m ~ r[e*p{- ((4”° d) ® - 4“ p, ) 2 /4 ■ (i) 

Here, x is the 14-dimensional vector describing a point 
in parameter space and z t are principal components of 
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the observables where each observable yi is first scaled 
by ai which describes the uncertainty one assigns to the 
comparison of the model to experiment, with cr,; account¬ 
ing for both experimental uncertainties and the error 
one might associate with the model missing some of the 
physics. Here, the uncertainties were all chosen to be 
6% of each observable. Changing this to 9% only mod¬ 
estly affected the final result. The largest source of un¬ 
certainty derives from the unknown impact of missing 
physics. These shortcomings will be discussed further 
below. 

Constraining the equation of state is the principal goal 
of this study. The equation of state was chosen to be 
consistent with that of a hadron gas for a tempera¬ 
ture of 165 MeV, which is the temperature at which the 
hydrodynamic description switched to the microscopic 
hadronic simulation. At the high-energy densities con¬ 
sidered here, one can neglect any small excess of baryons 
to antibaryons and the equation of state can be expressed 
in terms of a single variable such as the energy density 
e. For temperatures above 165 MeV, the speed of sound 
squared was parameterized to allow for a large range of 
equations of state, 


c s( e ) = c 2 A e h) + Q- 
X 0 = X'Rc s (e)Vl2, 


2/ A X 0 x + x2 

Cs[eh> ) X 0 x + x 2 +X'2 
x = In e/eh, 


( 2 ) 


where eh is the energy density corresponding to T = 165 
MeV. The two parameters R and X' describe the behav¬ 
ior of the speed of sound at energy densities above eh- 
Whereas R describes how the speed of sound rises or falls 
for small x, X' describes how quickly the speed of sound 
eventually approaches 1/3 at high temperature. Once 
given Cg(e), thermodynamic relations provide all other 
representations of the equation of state. Runs were per¬ 
formed for 0.5 < X' < 5, and with —0.9 < R < 2. In the 
limit R —> — 1 the speed of sound will have a minimum 
of zero. 

Ten of the 14 model parameters described the initial 
stress-energy tensor and flow used to describe the ini¬ 
tial state and instantiate the hydrodynamic calculation, 
with 5 separate parameters describing the initial state for 
each beam energy. Three parameters varied the trans¬ 
verse profile of the initial energy density at each beam 
energy: a weight between two saturation pictures, a nor¬ 
malization for the initial energy density, and a screening 
parameter. These three parameters, along with a param¬ 
eter used to vary the initial flow, are described in m- 
The fifth parameter describes the initial anisotropy of 
the stress energy tensor and was varied so that the lon¬ 
gitudinal pressure, T ZZl could vary between zero and the 
pressure P. The viscosity at T = 165 MeV and its tem¬ 
perature dependence were described by two parameters 
as was done in m and the final two parameters varied 
the equation of state. 

The details of both the physical model and the sta¬ 
tistical method are described in m The calculations 
shown here were based on 1200 full-model runs. Thirty 


observables, 15 for RHIC data and 15 for the LHC, were 
related to spectra, elliptic flow and femtoscopic source 
sizes. Observables were calculated for two centralities, 
20-30% and 0-5% for both the RHIC and LHC cases. At 
each centrality the spectral observables were the mean 
transverse momenta, (pt), for pions, kaons and protons, 
and the yield for pions. The three femtoscopic sizes, aver¬ 
aged over the experimentally analyzed momentum range, 
Rout, Rside and Riong described the dimensions of the 
outgoing phase space cloud of particles with the same 
momenta. The (pt)-weighted measurement of the elliptic 
flow, V '2 = (cos 2</>) quantified the preference for emitting 
particles in the reaction plane (<j> = 0 or 180°). Because 
the model used initial energy profiles that were smoothed 
by considering the averaged positions of incoming nucle¬ 
ons within a nucleus, rather than more realistic lumpy, 
or fluctuating, initial conditions, the model had to scale 
up its predictions for elliptic flow by a factor of 1.10. 
This accounts for the fact that the fluctuations result in 
larger initial transverse elliptic asymmetries which then 
lead proportionally to larger flows. The correction factor 
was quantitatively evaluated assuming a linear response 
in V 2 to initial eccentricity and found to be minimized 
in the 20-30% centrality class. The i >2 analysis was con¬ 
fined to the 20-30% centrality to minimize the effect of 
fluctuating initial conditions. 

The first 1000 runs were chosen semi-randomly 
throughout the 14-dimensional parameter space accord¬ 
ing to latin hyper-cube sampling. The thirty observables 
were then reduced to 14 principal components, which 
captured over 99.9% of the variance. Identically to what 
was done in D3!. these principal components were in¬ 
terpolated from the 1000 runs using a Gaussian process 
emulator during a Markov chain Monte Carlo (MCMC) 
exploration of the parameter space. This yielded a poste¬ 
rior sampling of the parameter space, i.e. a sampling that 
was weighted by the likelihood to reproduce the measured 
observables. A sampling of 50 points in parameter space 
was then chosen according to the posterior distribution 
and evaluated with the full model. Real model values 
were then compared to the emulated values at these 50 
points to validate the emulator in the regions of high 
likelihood which are most important in correctly deter¬ 
mining the posterior distribution. The emulator was then 
retrained using the 1050 runs and the validation proce¬ 
dure was repeated three additional times resulting in a 
total of 1200 model runs. The emulator’s accuracy in 
each case was found to be a few tenths of one unit when 
determining ln(£) in Eq. 0. The results shown here use 
emulation based upon the full model runs at these 1200 
points in parameter space, 200 of which are distributed 
according to the posterior distribution. 


III. RESULTS 

The ability of the procedure to accurately identify 
likely regions of parameter space is illustrated in Figs. 
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FIG. 1. Twenty pion, kaon and proton spectra as measured 
by the ALICE collaboration at the LHC (circles for 0-5% and 
squares for 20-30%) [14] are compared to model predictions 
using parameters randomly taken from the prior parameter 
space (panels a-c) and using parameters weighted by the like¬ 
lihood (d-f). 
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FIG. 2. The pion azimuthal anisotropy V 2 , often referred 
to as elliptic flow, from ALICE [16] for the 20-30% centrality 
(circles) are compared to model predictions using parame¬ 
ters randomly taken from the prior parameter space (a), and 
weighted by the likelihood (b). 


FIG. 3. Two-particle femtoscopic source sizes from ALICE 
[15] (circles for 0-5% and squares for 20-30% centrality) are 
compared to model predictions using parameters randomly 
taken from the prior parameter space (a-c), and weighted by 
the likelihood (d-f). 
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EE and [3] by comparing both full model calculations at 
20 random points in parameter space and then again at 
20 points chosen proportional to the likelihood defined in 
Eq. ( 0 . Calculations are compared to ALICE Collabo¬ 
ration at the LHC. Similarly good representations of the 
experimental data are found for RHIC data, with results 
very similar to those in m- The procedure readily iden¬ 
tified regions of parameter space that matched the exper¬ 
imental data within the 6% uncertainty assumed here. 
Nonetheless, it appears that the procedure finds spec¬ 
tra that have transverse momenta that are a few percent 
higher than the the experiment, and femtoscopic source 
sizes that are a few percent larger. This suggests the fem¬ 
toscopic data and the spectra are competing for agree- 


FIG. 4. The posterior likelihood for the two parameters that 
describe the equation of state, X' and R, have a preference 
to be along the diagonal. This shows that experiment con¬ 
strains some integrated measure of the overall stiffness of the 
equation of state, i.e. a softer equation of state just above T c 
is consistent with the data if it is combined with a more rapid 
stiffening at higher temperature. 

ment with the model, as slightly more explosive models 
would better match the femtoscopic observations, while 
less explosive models would better reproduce the spec¬ 
tra. This implies that improved physics might be needed 
if one were to reproduce the experimental results much 
better than 6%. 

The ability of the procedure to constrain the two pa- 
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FIG. 5. (a) Fifty equations of state were generated by ran¬ 

domly choosing X' and R in Eq. (|2| from the prior distribu¬ 
tion and weighted by the posterior likelihood (b). The two 
upper thick lines in each figure represent the range of lat¬ 
tice equations of state shown in [4], and the lower thick line 
shows the equation of state of a non-interacting hadron gas. 
This suggests that the matter created in heavy-ion collisions 
at RHIC and at the LHC has a pressure that is similar, or 
slightly softer, to that expected from equilibrated matter. 


rameters that determine the equation of state is shown 
in Fig. [4] As a function of X' and R defined in Eq. 
([ 2 ]), the likelihood is significant for a large band near 
the diagonal. Higher values of X ', which delays the ap¬ 
proach of the speed of sound to one third until higher 
energy densities and makes the equation of state softer, 
can be compensated by higher values of R, which sends 
the speed of sound higher just above T c and makes the 
equation of state stiffer. Fifty values of X' and R were 
then taken randomly from both the prior, and weighted 
by the posterior likelihood. For each case the speed of 
sound is plotted as a function of the temperature in Fig. 
[2] It is clear that the experimental results significantly 
constrain the equation of state and we also note that 
the RHIC and LHC data in combination provide a bet¬ 


ter constraint than either can alone. It appears that 
the speed of sound cannot fall much below the hadron 
gas value, ~ 0.15, for any extended range and that it 
must rise with temperature. Figure [5] also shows a range 
of equations of state from lattice calculations [3j. The 
equations of state found here show a preference for being 
slightly softer than those from the lattice, but the ranges 
overlap. 

IV. CONCLUSIONS 

Determining the equation of state from experiment 
has proven difficult due to the intertwined links between 
model parameters and numerous observables. The sta¬ 
tistical techniques applied here overcome these difficul¬ 
ties. The resulting constraints suggest the speed of sound 
gradually rises as a function of temperature from the 
hadron gas value. The band of equations of state from 
Fig. [5] is modestly softer than that of lattice calculations, 
but has significant overlap. This analysis strengthens the 
supposition that the matter created in relativistic heavy 
ion collisions has properties similar to that of equilibrated 
matter according to lattice calculations and shows that 
our model describes the dynamics of heavy ion collisions 
well enough to permit the extraction the thermodynamic 
and transport properties of equilibrium condensed QCD 
matter. 
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